clear all; close all; clc
 projectdir   = '/Volumes/NO NAME';
 projectdirDB = '/Users/tamoni/Dropbox';
 addpath(fullfile(projectdir, '/Backup/Papers/IRF Linear and NonLinear/Jorda/SLP/slpMatlab/LibrarySLP'));
%----------- Load shocks using two identification schemes ----------------%
LogVsLevel = 1; 
  NumberConsecutiveshocks = 2; % 1 previous shocks that is  positive, and the 2nd one is positive too 
% NumberConsecutiveshocks = 3; % 2 previous shocks that are positive, and the 3rd one is positive too 
% NumberConsecutiveshocks = 4; % 3 previous shocks that are positive, and the 4th one is positive too   

% Settings
%COUNTRY LABELS:
clabel ={'USA','Japan','Germany','Italy','UK','France','Canada','Spain',...
    'Sweden','Switzerland','Netherlands'};

%WHICH MACRO UNCERTAINTY INDEX? Baker, Bloom & Davis (2016) (BBD) OR JLN
muindex=0; % 1 =Fin U ; 0 = BBD

% G7 + Spain + Sweden + Switzerland + Netherlands
 % 1.USA, 2. JP , 3. GER, 
 % 4. IT, 5. UK,  6. FRA, 
 % 7. CN, 8. ESP, 9. SWE, 10. CHE, 11. NL
country = 5; %2,...,11

LoadShocksMandQpositiveEPUrevisedMAC; %sample=[ 1961 1; 2019  12];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%DATA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
data = dataMatrix;
% % startsample=find(data(:,1)==sample(1,1)*100+sample(1,2));
% % endsample  =find(data(:,1)==sample(2,1)*100+sample(2,2));

% % data = data(startsample:endsample,2:end);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%DATA TRANSFORMATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % data(:,[3:10]) = log(data(:,[3:10]))*100;
%--------------------------------------%
clear startsample endsample 

T = size(data,2);
 P  =  2; %                including in the LP regression P lags of all variables in the system
%P  =  4; % 
 
% plottype = 'riskfr';      ii=3;                            
  plottype = 'Output';      ii=3; %
%                           ii=4; % for hours worked
% plottype = 'Inflat';     %ii=?; % 
% plottype = 'Credit' ;    %ii=2; %      
%% IRF 
h1 = 0 ;  %start LP at h=0 or h=1 (h=1 if impose no contemporaneous impact)
H  = 4*3; %#horizon of IR

% We are interested in the estimation of the dynamic multiplier of y_t+h with respect to a change in x_t for h ranging from 0 to H
%------------ IRF of GDP growth to a shock variable ----------------------%
 y  = data(:,ii);    %endogenous variable
% if ii == 8 || ii == 9
%     y  = [NaN(12,1); data(13:end,ii)-data(1:end-12,ii)];    %endogenous variable
% end
%-------------------------------------------------------------------------%
% Remove_Large_Shocks = 1;
    x  = uEPU;      
%/////////////////////////////////////////////////////////////////////////%
% w  = [lagmatrix( data(:,unique([ii 3       ])  ), 1:P )                       ]; %control variables (contemporaneous vars, lagged vars)
% w  = [lagmatrix( data(:,unique([ii 4       ])  ), 1:P )                       ]; %control variables (contemporaneous vars, lagged vars)
% w  = [lagmatrix( data(:,unique([ii 2       ])  ), 1:P ) lagmatrix(logEPU_rescaled,1:P )];
  w  = [lagmatrix( data(:,unique([ii 2  4    ])  ), 1:P ) lagmatrix(logEPU_rescaled,1:P )];
OKforReg = isfinite(x) & isfinite(y);
for cc=1:size(w,2)
    OKforReg = OKforReg & isfinite(w(:,cc)) ;
end  
  x = x(OKforReg,:); 
  y = y(OKforReg,:)*100; 
  if country == 7 || country == 8 || country == 11
      y = y/3; 
  end
  w = w(OKforReg,:);
  
% This shrinks the estimated dynamic multiplier to a polynomial of order r - 1
r = 3; %(r-1)=order of the limit polynomial (so r=3 implies the IR is shrunk towards a quadratic polynomial)

%-------------------------------------------------------------------------%
% k-fold cross validation 
  k_fold=5;
% slp = locproj(y,x,w,h1,H,'smooth',r,0.01);
% 
% lambda = [ 1:0.5:10] * T;
% slp    = locproj_cv(slp,k_fold,lambda);
% 
% figure(2)
% plot( lambda , slp.rss , '-o' )
% 
% lambda_opt = lambda( min( slp.rss ) == slp.rss );
% 
% fprintf('Lambda from 5-fold cross validation is %3.3f \n',lambda_opt);
% %lambda = lambda_opt;
%-------------------------------------------------------------------------%       
                 lambda =  20; %value of penalty 
if country == 5; lambda =  2000; end %value of penalty
 lp = locproj(y,x,w,h1,H,'reg');             %IR from (standard) Local Projection
slp = locproj(y,x,w,h1,H,'smooth',r,lambda); %IR from Smooth Local Projection

figure(41)
%%
 scalingPLOTS = 1.0*volaShock;
%scalingPLOTS = 3.0; %Remember that BBD compute model-implied responses of industrial production and employment to a 90-point upward EPU innovation
  plot( 0:H,scalingPLOTS*[ lp.IR slp.IR],'LineWidth',3); hold on;  
% plot( 0:H,scalingPLOTS*100*  lp.Interval_up,'r--');           plot( 0:H ,scalingPLOTS* lp.Interval_down,'r--');
  plot( 0:H,scalingPLOTS* slp.Interval_up,'r--','LineWidth',2); plot( 0:H ,scalingPLOTS*slp.Interval_down,'r--','LineWidth',2);
plot( 0:H , zeros(H+1,1) , '-k' , 'LineWidth' , 2 )
grid
xlim([0 H])

switch plottype
    case 'Output' 
        title('Output','FontSize',16);     %set(gca,'XTick',t(4:4:end),'FontSize',12);
       %ylim([-1.6 1.0]);set(gca,'YTick',[-1.6:0.5:1.6],'FontSize',12);ylabel('Percent','FontSize',12)
%     case 'riskfr'
%         title('Risk-free rate','FontSize',16);%set(gca,'XTick',t(4:4:end),'FontSize',12);
%         %ylim([-0.4 0.4]);set(gca,'YTick',[-0.6:0.2:0.6],'FontSize',12);ylabel('Percent','FontSize',12)
%     case 'SP500'
%         title('Stock market','FontSize',16); %set(gca,'XTick',t(4:4:end),'FontSize',12);
%         ylim([-4.0 1.4]);set(gca,'YTick',[-4.0:0.5:1.4],'FontSize',12);ylabel('Percent','FontSize',12)
%     case 'Inflat'
%         title('Inflation','FontSize',16); %set(gca,'XTick',t(4:4:end),'FontSize',12);
%         %ylim([-5.4 5.4]);set(gca,'YTick',[-5.2:0.5:5.2],'FontSize',12);ylabel('Percent','FontSize',12)    
%     otherwise
        warning('Unexpected plot type. No plot created.')
end
legend('IR_{lp}','IR_{slp}',                             'Location','Best')

keep slp ii 
fname = sprintf('EPUBBD_LP_linear.mat');
%if     ii==5
    slpOutput     = slp.IR; slpOutputUP     = slp.Interval_up; slpOutputDown     = slp.Interval_down; 
    save(fname,'slpOutput','slpOutputUP','slpOutputDown')
%end    